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SUMMARY 

We  extend  the  explicit  in  time  high-order  triangular  discontinuous  Galerkin  (DG)  method  to  semi- 
implicit  and  then  apply  the  algorithm  to  the  two-dimensional  oceanic  shallow  water  equations;  we 
implement  high-order  semi-implicit  time-integrators  using  the  backward  difference  formulas  from 
orders  one  through  six.  The  reason  for  changing  the  time-integration  method  from  explicit  to  semi- 
implicit  is  that  the  explicit  method  requires  a  very  small  time-step  in  order  to  maintain  stability, 
especially  for  high-order  DG  methods.  Changing  the  time-integration  method  to  be  semi-implicit 
allows  one  to  circumvent  the  stability  criterion  due  to  the  gravity  waves,  which  for  most  shallow  water 
applications  are  the  fastest  waves  in  the  system  (the  exception  being  supercritical  flow  where  the 
Froude  number  is  greater  than  one).  The  challenge  of  constructing  a  semi-implicit  method  for  a  DG 
model  is  that  the  DG  machinery  requires  not  only  the  standard  finite  element-type  (FE)  area  integrals 
but  finite  volume-type  (FV)  boundary  integrals  as  well.  These  boundary  integrals  pose  the  biggest 
challenge  in  a  semi-implicit  discretization  because  they  require  the  construction  of  a  Riemann  solver 
that  is  the  true  linear  representation  of  the  nonlinear  Riemann  problem;  if  this  condition  is  not  satisfied 
then  the  resulting  numerical  method  will  not  be  consistent  with  the  continuous  equations.  In  this  paper 
we  present  semi-implicit  time- integrators  for  the  DG  method  that  maintain  most  of  the  usual  attributes 
associated  with  DG  methods  such  as:  high-order  accuracy  (in  both  space  and  time),  parallel  efficiency, 
excellent  stability,  and  conservation.  The  only  property  lost  is  that  of  a  compact  communication  stencil 
typical  of  time-explicit  DG  methods;  implicit  methods  will  always  require  a  much  larger  communication 
stencil.  We  apply  the  new  high-order  semi-implicit  discontinuous  Galerkin  method  to  the  shallow  water 
equations  and  show  results  for  many  standard  test  cases  of  oceanic  interest  such  as:  standing,  Kelvin 
and  Rossby  soliton  waves,  and  the  Stommel  problem.  The  results  show  that  the  new  high-order  semi- 
implicit  DG  model,  that  has  already  been  shown  to  yield  exponentially  convergent  solutions  for  smooth 
problems,  results  in  a  more  efficient  model  than  its  explicit  counterpart.  Furthermore,  the  capacity 
to  use  high-order  time-integrators  offers  a  big  advantage  in  accuracy  when  simulating  time-dependent 
problems  especially  when  using  high-order  DG  methods;  without  high-order  time-integration  it  makes 
little  sense  to  use  high-order  spatial  discretizations.  Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 
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1.  Introduction 

The  discontinuous  Galerkin  (DG)  method  has  come  into  prominence  in  the  last  decade  in  all 
areas  of  numerical  modeling;  however,  it  has  only  been  in  the  last  few  years  that  this  method  has 
received  attention  in  geophysical  fluid  dynamics.  The  high-order  accuracy,  geometric  flexibility 
to  use  unstructured  grids,  conservation,  and  monotonicity  properties  of  the  DG  method  makes 
it  a  prime  candidate  for  the  construction  of  future  ocean  and  shallow  water  models.  The 
advantages  offered  by  the  DG  method  will  benefit  all  areas  of  ocean  modeling  but,  specifically, 
it  will  improve  coastal  ocean  models  where  proper  coastline  representation,  and  the  ability  to 
handle  steep  gradients  should  translate  into  better  modeling  of  tsunamis,  storm  surges,  and 
hurricanes.  Let  us  now  review  the  literature  concerning  the  application  of  the  DG  method  to 
the  oceanic  shallow  water  equations. 

Schwaneberg  and  Kongeter  (2000)  [32]  first  used  the  DG  method  for  the  planar  shallow  water 
equations,  followed  by  the  work  of  Li  and  Liu  (2001)  [26],  and  Aizinger  and  Dawson  (2002) 
[2],  Dupont  and  Lin  (2004)  [12],  Eskilsson  and  Sherwin  (2004)  [13],  Remade  et  al.  (2006)  [28], 
and  Kubatko  et  al.  (2006)  [25]  constructed  shallow  water  models  on  triangles  using  a  collapsed 
local  coordinate  discontinuous  Galerkin  method.  Giraldo  and  Warburton  (2008)  [20]  developed 
a  high-order  DG  oceanic  shallow  water  model  on  unstructured  adaptive  triangular  grids.  In 
that  paper,  we  showed  that  the  model  yields  exponentially  convergent  solutions  (for  smooth 
problems).  However,  we  used  explicit  time- integration  methods  which,  while  easy  to  implement, 
require  small  time-steps  in  order  to  maintain  stability.  To  ameliorate  this  deficiency  found  in 
all  DG  shallow  water  models,  we  extend  the  explicit  time-integrators  to  semi- implicit.  To 
date,  there  has  been  no  work  on  the  development  of  semi-implicit  time-integrators  for  shallow 
water  DG  models;  all  of  the  DG  shallow  water  models  found  in  the  literature  use  explicit 
time-stepping,  including  those  discussed  above.  Furthermore,  the  only  work  on  DG  and  semi- 
implicit  methods  found  in  the  literature  are  the  papers  by  Dolejsi,  Feistauer,  and  co-authors  on 
the  compressible  Navier-Stokes  equations  (see  [11],  [10],  [14],  [15],  [9]).  Their  semi- implicit  DG 
formulation  is  based  on  low-order  polynomial  spaces  (third  order  or  less)  and  their  approach  is 
fundamentally  different  from  ours  in  that  they  rely  on  a  linearization  of  the  nonlinear  operators 
in  conjunction  with  a  special  flux  function  that  facilitates  this  linearization.  Our  approach  [29] 
relies  on  extracting  the  linear  operators  containing  the  fastest  wave  speeds  in  the  system 
and  then  discretizing  them  implicitly  in  time.  While  both  approaches  are  very  effective,  our 
approach  is  more  similar  to  the  classical  semi-implicit  method  first  proposed  by  Robert  et  al. 
(1972)  [31].  The  advantage  of  this  approach  is  that,  once  the  numerical  machinery  is  developed, 
it  can  be  applied  to  any  equation  set  with  minimal  modifications.  Moreover,  the  semi-implicit 
DG  approach  is  easily  extendable  to  generalized  families  of  linear  multi-step  time-integration 
methods  as  we  show  here. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  describes  the  governing 
equations  of  motion  used  to  test  our  numerical  method.  In  Sec.  3  we  describe  the  spatial 
discretization  of  the  governing  equations  and  in  Sec.  4  the  time-integrators  used.  Finally,  in 
Sec.  5  we  present  comparisons  between  the  explicit  and  semi-implicit  versions  of  the  model. 
This  then  leads  to  a  summary  on  the  direction  of  future  work. 
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2.  Continuous  Equations 


The  oceanic  shallow  water  equations  are  a  system  of  nonlinear  partial  differential  equations 
which  govern  the  motion  of  a  viscous  incompressible  fluid  in  a  shallow  depth.  The  predominant 
feature  of  this  type  of  fluid  is  that  the  characteristic  length  of  the  fluid  is  far  greater  than  its 
depth,  this  is  analogous  to  ocean  flow  problems  and  is  the  reason  these  equations  are  typically 
used  as  a  first  step  toward  the  construction  of  ocean  models. 

The  shallow  water  equations  in  conservation  form  are 

+  V  •  F(q)  =  S(q)  (1) 

where  q  =  (<j),  UT)T  are  the  conservation  variables, 


is  the  flux  tensor  and 

S{q)  =  ~  (  .f(kxU)-  ct>sVct>B 

is  the  source  function  where  the  nabla  operator  is  defined  as  V  =  ( dx,dy)T ,  ®  denotes  the 
tensor  product  operator,  <j)  =  g(hs  +  hs )  is  the  geopotential  height  where  g  is  the  gravitational 
constant,  hs  is  the  free  surface  height  of  the  fluid,  (j>B  =  ghs  is  the  bathymetry  (e.g.,  bottom 
of  the  ocean)  which  we  assume  constant,  U  =  <f>u  is  the  momentum,  u  =  (u,  v)T  is  the  velocity 
vector,  /  =  /o  +  (3(y  —  ym )  is  the  Coriolis  parameter,  k  =  (0,  0, 1)T  is  the  unit  normal  vector 
of  the  x-y  plane,  and  the  term  X-2  is  a  rank- 2  identity  matrix.  The  vector  t  is  the  wind  stress, 
p  is  the  density,  H  is  a  mean  height,  and  the  constant  7  is  the  bottom  friction. 


2.1.  Linearized  Continuous  Equations 


Let  us  now  decompose  Eqs.  (1)  -  (3)  into  their  linear  and  nonlinear  components.  Splitting  the 
geopotential  height  (j>  into  the  depth  from  mean  sea  level  to  the  ocean  bottom  <f >b  and  the 
height  from  mean  sea  level  to  the  water  surface  4>s  we  then  have  <j>(x,t)  =  <f>s(x,t)  +  4>b(x) 
which  we  can  now  use  to  substitute  into  the  equations  to  get 


d0s 

dt 


dU 

dt 


V  • 


Snl 


U®U 


+  2  )  +  4>s4>b^2 


V  •  U  =  0  (4) 

=  -/(fcx[/)  +  <feV0B  (5) 


pH 


where  5nl  is  a  switch  which  retains  the  nonlinear  terms  when  5nl  =  1  and  turns  them  off  for 
Snl  =  0. 

In  the  next  few  sections,  we  describe  the  semi-implicit  (SI)  method  for  the  oceanic  shallow 
water  equations  in  the  particular  case  that  the  discontinuous  Galerkin  method  is  used  to 
represent  spatial  derivatives.  As  the  reader  will  see,  a  pivotal  component  of  the  SI  method  is 
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the  construction  of  a  linearized  form  of  the  continuous  equations.  Linearizing  Eqs.  (1)  -  (3) 
yields 


d(j)s 

dt 


V-U  =  0 


dU 

+  V  •  ((/>s<^bX2) 


-/(fe  x  C7)  +  +  —  -  7C/ 

pH 


(6) 


which  are  obtained  by  setting  Snl  =  0  in  Eqs.  (4)  and  (5).  The  maximum  eigenvalue  of  the 
linear  system  given  in  Eq.  (6)  is  A l  =  \Z4>b  which  is  in  fact  the  linearized  eigenvalue  of  A a tl 
obtained  for  the  nonlinear  system  given  in  Eqs.  (1)  -  (3).  From  Eq.  (6)  we  can  define  the 
linear  operator  as  follows 


L  =  _(  V-£7 

V  V  •  (^fel2)  +  f(k  xU)-  </>s  V 4>b  +  7 u 

We  will  return  to  this  linear  operator  in  Sec.  4.  Let  us  now  describe  the  approximation  of 
the  spatial  derivatives  by  the  DG  method.  We  need  to  know  this  information  before  we  can 
construct  the  semi-implicit  solution. 


3.  Triangular  Discontinuous  Galerkin  Method 

In  this  section  we  describe  the  approximation  of  the  spatial  derivatives  of  the  shallow  water 
equations  using  the  discontinuous  Galerkin  method  on  triangles.  This  includes:  the  choice  of 
basis  functions,  integration,  construction  of  the  semi-discrete  problem  and  its  corresponding 
matrix  form. 

3.1.  Basis  Functions 

To  define  the  discrete  local  operators  we  begin  by  decomposing  the  domain  O  into  Ne 
conforming  non-overlapping  triangular  elements  f le  such  that 

Ne 

n=\Jne. 

e=l 

The  condition  on  grid  conformity,  however,  is  not  required  by  the  DG  method;  we  only  impose 
this  condition  to  simplify  the  discussion. 

To  perform  differentiation  and  integration  operations,  we  introduce  the  nonsingular  mapping 
x  =  which  defines  a  transformation  from  the  physical  Cartesian  coordinate  system 

x  =  ( x,y)T  to  the  local  reference  coordinate  system  £  =  ( £,ri)T  defined  on  the  reference 
triangle  fie  =  {(£,77),  -1  <  £,r?  <  1,  £  +  77  <  0,}. 

Let  us  now  represent  the  local  element-wise  solution  q  by  an  JVth  order  polynomial  in  £  as 

Mn 

<1n(€)  =  '52i,i(€)QN(€i)  (8) 

i=  1 

where  represents  Mjv  =  ^(N  +  1){N  +  2)  interpolation  points  and  '0i(£)  are  the  associated 
multivariate  Lagrange  polynomials.  For  the  interpolation  points  we  choose  the  nodal  sets 
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based  on  the  electrostatics  [21]  and  Fekete  [36]  points;  for  simplicity  we  shall  refer  to  these 
nodal  sets  collectively  as  Fekete  points.  We  have  already  described  the  construction  of  the 
nodal  basis  functions  in  [18],  [20]  and,  for  the  sake  of  brevity,  omit  this  discussion  here. 


3.2.  Integration 


3.2.1.  Area  Integrals  In  order  to  complete  the  discussion  of  the  local  element-wise  operations 
required  to  construct  discrete  spatial  operators  we  must  describe  the  integration  procedure 
required  by  the  integral  formulation  of  all  Galerkin  methods.  For  any  two  functions  /  and  g 
the  2D  (area)  integration  I  a  proceeds  as  follows 


Wg] 


where  Mq  is  a  function  of  Q  which  represents  the  order  of  the  cubature  approximation.  For 
Wi  and  we  use  the  high-order  cubature  rules  for  the  triangle  given  in  [35,  6,  27,  7];  because 
we  use  order  2N  integration,  which  is  exact  for  this  equation  set,  then  neither  spatial  filters 
nor  smoothing  diffusion  operators  are  used.  Furthermore,  we  omit  the  use  of  slope  limiters. 


3.2.2.  Boundary  Integrals  The  DG  method  also  requires  the  evaluation  of  boundary  integrals 
which  is  the  mechanism  by  which  the  fluxes  across  element  edges  are  evaluated  and  allows  the 
discontinuous  elements  to  communicate.  For  any  two  functions  /  and  g  the  ID  (boundary) 
integration  2 b  proceeds  as  follows 

/'  Q 

ZB[f,g\  =  /  f(x)g(x)dx  =  I  Js(£,;)  I  f(£i)g(£i) 

*= o 

where  Q  represents  the  order  of  the  quadrature  approximation.  Using  Legendre-Gauss 
quadrature  we  can  use  Q  =  N  to  achieve  order  2 N  accuracy. 


3.3.  Tangent  and  Normal  Vectors  of  the  Element  Edges 

Below  it  will  become  evident  that  in  order  to  construct  a  discontinuous  Galerkin  discretization 
requires  knowing  a  bit  about  the  element  geometry.  In  continuous  Galerkin  methods,  such  as 
the  finite  element  method,  the  only  required  information  is  the  basis  functions,  metric  terms, 
and  cubature  rules.  The  DG  method  requires  all  of  this  finite  element-type  information  plus 
some  finite  volume-type  information  regarding  the  element  edges  and  the  element  neighbors 
sharing  these  edges.  However,  the  good  news  for  the  DG  method  is  that  regardless  of  the  order 
of  the  basis  function,  N ,  each  element  only  has  three  edge  neighbors  (this  is  true  only  for 
conforming  grids).  This  is  the  process  by  which  a  DG  element  shares  its  local  information  with 
its  neighbors. 


3.f.  Semi-Discrete  Equations 

Applying  the  discontinuous  Galerkin  discretization  to  Eq.  (1),  and  using  Green’s  theorem 
yields  the  classical  DG  which  we  refer  to  as  the  weak  form 

f  (  ~dt  ~  '  V  _  '5^)N)  ^i(x)dx  =  ~  /  ^(*)  n(e’°  ‘  FNl)dx  (9) 

J \  J  1=1 
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where  Fjv  =  F(qN)  and  Sn  =  S(qN)  with  F  and  S  given  by  Eqs.  (2)  and  (3),  respectively. 
Note  that  Eq.  (9)  states  that  qN  satisfies  the  equation  on  each  element  fle  for  all  if)  £  S  where 
S  is  the  finite-dimensional  space 

S  =  {tp  €  C2^F)  :  '0|ne  €  Pjv(f2e)Vf2e}  , 


PN  is  the  polynomial  space  defined  on  f \e  and  the  union  of  these  elements  defines  the  entire 
global  domain  -  that  is,  f \  =  U^i  with  Ne  representing  the  total  number  of  triangular 
elements.  It  should  be  mentioned  that  in  DG  methods,  the  space  Pn~Pn  can  be  used  without 
violating  the  inf-sup  (Ladyzlrenskaya-Babuska-Brezzi)  condition  which  must  be  observed  by 
continuous  Galerkin  methods  (such  as  the  finite  element  method)  in  order  to  avoid  the  effects 
of  spurious  pressure  modes. 

In  the  boundary  integral  of  Eq.  (9)  n  is  the  outward  pointing  unit  normal  vector  of  the 
element  edge  Te  and  F*N  is  the  Rusanov  numerical  flux  (although  other  fluxes  are  also  possible) 


p(*’1) 

r  N 


+  fn  (<$)  -|a«i(4} 


(10) 


where  A®  =  max  +  \/ 4>^e\  \U^\  +  \J with  t/(e>0  =  u^e’1^  ■  n W  being  the  normal 

component  of  velocity  with  respect  to  the  edge  Te,  and  the  superscripts  e  and  l  represent 
the  element  e  and  its  three  edge  neighbors  l.  The  normal  vector  is  defined  as  pointing 

outward  from  the  element  e  to  its  edge  neighbor  l. 

Integrating  Eq.  (9)  by  parts  once  more  yields  the  strong  form 


( 

^  dt 


dx 


(11) 


which,  although  mathematically  equivalent  to  the  weak  form,  yields  different  numerical 
solutions.  Based  on  previous  studies  (see  Giraldo  [18]  and  Kopriva  [24])  we  use  the  strong 
form  exclusively  in  this  paper. 


3.5.  Matrix  Form  of  the  Semi-Discrete  Equations 
Substituting  the  polynomial  approximation 


Mjv 

»= i 

into  Eq.  (11)  we  can  now  write  the  semi-discrete  system  as 


ipiipj  dx  dx  -Ff-  f  dx  Sje)  =  Y,  [  ^ j"{e,,)  dx  ■  (F(6)  -  pi*'l) 

rig  rig  | ^  r g 

(12) 

Next,  note  that  by  defining  the  following  element  matrices 

Mif  =  I  4>if)jdx,  Mffl)  =  f  ipiipjn(-e’l'>  dx,  D*f)  =  f  tfiVifjdx 

« rig  «/  r  g  « rig 
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allows  us  to  write  Eq.  (12)  in  the  following  matrix  form 


(e) 


M<e>— L. 
'3  dt 


- 


E  (Mif)T  (F(e) 


(13) 


where  the  superscript  e  denotes  an  element-based  evaluation  and  l  denotes  edge-based 
evaluations.  Next,  using  the  approach  described  in  [20]  for  eliminating  the  mass  matrix,  we 
write 


(e) 

d  = 


(M(e))  1  D{e\  M{e'l)  = 
which  then  allows  us  to  write  Eq.  (13)  as  follows 


dq 


(e) 


dt 


(A?)  Ff  -  Sie)  =  E  (MS’°)  (f(e)-FM)/  (!4) 


Equation  (14)  is  the  form  we  shall  use  in  the  construction  of  the  explicit  and  semi-implicit 
discretizations. 


3.6.  Boundary  Conditions 

In  all  the  test  cases  we  only  consider  no-flux  boundary  conditions;  we  will  extend  our  model 
to  more  general  boundary  conditions  in  future  work.  The  no-flux  boundary  conditions  are 
enforced  by  virtue  of  the  statement 

n  ■  u  =  0  (15) 

at  the  boundaries.  Thus,  we  seek  to  eliminate  the  normal  component  of  the  velocity  to  the  no¬ 
flux  boundary  without  altering  the  tangential  component.  The  tangent  vector  to  a  boundary 
is  obtained  by  t  =  k  x  n  which  is  equal  to  t  =  — nyi  +  nxj.  Thus  we  solve  the  following  2x2 
system: 


where  ut  =  t  ■  u  is  the  tangential  component  of  velocity.  This  boundary  condition  is  imposed 
only  weakly  through  the  boundary  integrals  in  Eq.  (14);  that  is,  it  only  comes  in  through  the 
Rusanov  flux. 


4.  Time-Integrator 

In  Sec.  3  we  described  the  approximation  of  the  spatial  derivatives  using  the  DG  method.  We 
are  now  in  a  position  to  describe  the  approximation  of  the  time  derivatives.  Let  us  begin  with 
the  description  of  the  explicit  time-integration  followed  by  the  semi-implicit  time-integration 
methods.  We  use  a  method  of  lines  approach  for  both  methods. 

4-1-  Explicit  Method 

In  order  to  advance  the  solution  in  time  while  retaining  some  level  of  high-order  accuracy 
we  use  third  order  strong  stability  preserving  (SSP)  Runge-Kutta  methods  (see  [5]  and  [33]). 
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For  completeness  we  define  them  now.  Let  us  write  the  semi-discrete  (in  space)  equations  as 
follows 

!-*«> 

where 


i(«) 


T 


R(q)  =  -lD£'  F(ie)(q)  +  Sr  + 


?(e) 


3 

E 

1=1 


M 


(e,0 


F(e)(q)  - 


(*.0/ 


(17) 


and  q  is  in  fact  q^e\  that  is,  the  solution  within  each  element  e. 
The  SSP  temporal  discretization  of  this  semi-discrete  equation  is 


for  k  =  1, S  : 

qk  =  ak0qn  Fa^-1 +ak2qk~3 +f3kAtR(qk-1) 

where  q°  =  qn,  qs  =  qn+1,  S  denotes  the  number  of  RK  stages,  and  the  coefficients  a  and  /? 
are  given  in  Table  I. 


Method 

k 

ot  0 

Ol2 

P 

RK3 

1 

1 

0 

0 

1 

2 

3/4 

1/4 

0 

1/4 

3 

1/3 

2/3 

0 

2/3 

RK34 

1 

1 

0 

0 

1 

2 

0 

1 

0 

1/2 

3 

2/3 

1/3 

0 

1/6 

4 

0 

1 

0 

1/2 

RK35 

1 

1 

0 

0 

0.377268915331368 

2 

0 

1 

0 

0.377268915331368 

3 

0.355909775063327 

0.644090224936674 

0 

0.242995220537396 

4 

0.367933791638137 

0.632066208361863 

0 

0.238458932846290 

5 

0 

0.762406163401431 

0.237593836598569 

0.287632146308408 

Table  I.  Coefficients  for  the  explicit  strong  stability  preserving  third  order  Runge-Kutta  methods. 


In  Fig.  1,  we  show  the  stability  regions  of  these  methods  for  the  equation 


dq 

dt 


=  A  q 


which  is  a  proxy  for  an  advection-dominated  equation;  a  reasonable  approximation  for  the 
shallow  water  equations.  Figure  1  shows  that  the  stability  region  of  RK35  is  larger  than  those 
for  RK34  and  RK3.  We  are  interested  specifically  in  the  region  along  the  imaginary  axis 
because  this  is  the  region  most  important  to  advection-dominated  problems.  In  [8]  all  three 
RK  methods  were  studied  for  the  Navier-Stokes  equations  and  it  was  determined  that  RK35 
is  indeed  the  most  efficient  of  the  third  order  methods.  We  have  found  similar  results  here 
for  the  shallow  water  equations.  Based  on  the  results  of  these  studies,  we  use  RK35  for  the 
comparisons  with  the  semi-implicit  methods  which  we  now  describe. 
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Re(z) 


Figure  1.  The  stability  regions  for  the  explicit  strongly  stability  preserving  third  order  Runge-Kutta 
methods,  RK3  (three  stage),  RK34  (four  stage),  and  RK35  (five  stage). 


4-2.  Semi- Implicit  Method 

To  extend  the  size  of  the  time-step  we  use  a  generalized  semi-implicit  method  of  order  K .  Let 
us  write  Eqs.  (l)-(3)  in  the  following  compact  vector  form 

I  =  *<«>  <18> 

where  q  =  (<j)s,  UT)T  and  R(q)  is  defined  in  Eq.  (17).  This  system  can  be  represented  by  the 
equivalent  form 

^  =  {R(q)-6SIL(q)}  +  [6siL(q)}  (19) 

where  L(q)  is  the  linear  approximation  to  R  given  in  Eq.  (7)  and  contains  the  gravity  wave 
terms  (i.e. ,  the  fastest  waves  in  this  system,  at  least  for  subcritical  flow).  In  Eq.  (19)  the 
curly  brackets  denote  explicit  time-integration  while  the  square  brackets  denote  implicit  time- 
integration.  Note  that  the  variable  Ssi  is  a  switch  that  yields  a  fully  explicit  method  for  6si  =  0 
and  the  semi- implicit  method  for  5 si  =  1. 

As  was  done  in  [17,  19]  we  now  write  the  time-discretization  in  the  general  form: 

K  K 

qn+1  =  Y,  akqn~k  +  7  A  tYl3k[R(qn-k)  -  6SiL(qn~k )]  +  7  A  tSSiL(qn+1)  (20) 

k= 0  k= 0 

where  K  denotes  the  order  of  the  time-integrator.  To  simplify  the  discussion  of  the  semi-implicit 
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formulation,  let  us  now  introduce  the  following  auxiliary  variables 

K 

qtt  =  (21) 

k=0 

K  K 

qE  =  ^afcgn-fc+7At^/3fci?(qre-fc),  (22) 

k= 0  k=0 

K 

q  =  qE-J2a^n~k  (23) 

k= 0 

which  then  allows  us  to  write  Eq.  (20)  as 

Qtt=^  +  XL(Qtt)  (24) 

where  A  =  7 A tSsi-  In  Table  II  we  list  the  coefficients  for  the  backward  difference  formulas  of 
order  K  =  1, 6.  In  Figs.  2a  and  2b  we  show  the  stability  regions  of  the  explicit  and  implicit 


K= 

=  1 

K=2 

K=3 

K=4 

K=5 

K=6 

cto 

1 

4/3 

18/11 

48/25 

300/137 

360/147 

oil 

0 

-1/3 

-9/11 

-36/25 

-300/137 

-450/147 

Ct2 

0 

0 

2/11 

16/25 

200/137 

400/147 

Ct3 

0 

0 

0 

-3/25 

-75/137 

-225/147 

OL\ 

0 

0 

0 

0 

12/137 

72/147 

a  5 

0 

0 

0 

0 

0 

-10/147 

7 

1 

2/3 

6/11 

12/25 

60/137 

60/147 

A) 

1 

2 

3 

4 

5 

6 

01 

0 

-1 

-3 

-6 

-10 

-15 

02 

0 

0 

1 

4 

10 

20 

03 

0 

0 

0 

-1 

-5 

-15 

04 

0 

0 

0 

0 

1 

6 

05 

0 

0 

0 

0 

0 

-1 

Table  II.  Coefficients  for  the  backward  difference  formulas  of  orders  K  =  1, ...,  6. 

BDF  methods.  In  Fig.  2a  the  closed  loops  are  the  stability  regions  of  the  explicit  BDF  methods 
while  in  Fig.  2b  the  closed  loops  are  the  regions  of  instability  of  the  implicit  BDF  methods. 
For  the  shallow  water  equations,  we  are  interested  in  the  region  near  the  imaginary  axis  (for 
Re{z)  =  0). 

Let  us  now  describe  the  semi-implicit  method  in  terms  of  the  governing  equations.  Note 
that  the  operator  R  referenced  above  is  the  same  operator  described  for  the  explicit  method. 
However,  let  us  now  write  the  full  expression  given  in  Eq.  (24)  in  terms  of  the  operator  L. 
Substituting  the  linear  operator  defined  in  Eq.  (7)  into  Eq.  (24)  results  in  the  following  system 

<j>tt  =  <j>-  XV  -  Utt 

Utt  =  U-  A  [V  •  (<^B:z:2)  + /(fc  x  Utt) -fa Vc/>B+jUtt\  (25) 

where  we  have  retained  the  continuous  spatial  operators  for  clarity.  At  this  point,  there  is  no 
difference  between  the  semi-implicit  formulation  of  a  continuous  Galerkin  (e.g.,  finite  elements) 
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a)  b) 

Figure  2.  The  stability  regions  of  the  a)  explicit  and  b)  implicit  backward  difference  formulas  of  order 

K  =  1, ...,  6. 


and  a  discontinuous  Galerkin  method.  The  differences  arise  through  the  method  selected  for 
the  discretization  of  the  spatial  operators.  Replacing  the  continuous  spatial  operators  with  the 
DG  discrete  representations  yields 

r  3 

-  E 

1=1 


D 


(e) 


u 


(e) 

tt 


U[e)  -A  ID 


de) 


T 


(e) 


M 


Z=1 


Re.J) 


T 


(< ht(t>B^2 ) 


(MbZ 2)(e)  - 


'f(kxU%))-<f$V4>B+'yU$ 


(27) 


where  the  flux  values  are  specifically  defined  as 


U 


(*>0 

tt 


C7. 


(0 


(28) 


and 


+  {MbT^  -  | Ai| n(i)  (t/t(P  -  «7ite)) 


(29) 


These  equations  can  now  be  solved  as  a  coupled  system  of  linear  equations  for  </>tt  and  Utt. 
We  use  a  GMRES  solver  with  Jacobi  preconditioning  to  solve  this  system.  While  this  choice  of 
preconditioner  is  not  optimal,  the  resulting  iterative  method  is  nonetheless  robust  and  efficient 
in  terms  of  computational  time  and  memory  requirements  since  no  global  matrix  problem  ever 
needs  to  be  stored. 

In  standard  semi-implicit  methods  (e.g.,  see  [16],  [17],  [19]),  upon  writing  the  fully  discrete 
system  the  goal  is  then  to  apply  a  block  LU  decomposition  in  order  to  reduce  the  vector 
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system  of  equations  into  an  equivalent  scalar  system  of  equations;  for  first  order  systems  of 
equations,  the  resulting  problem  is  in  fact  quite  similar  to  a  Helmholtz  equation.  Constructing 
the  Helmholtz  problem  for  general  DG  polynomial  spaces  and  boundary  conditions  remains  an 
open  problem.  Thus  far,  we  only  know  how  to  construct  the  Helmholtz  problem  for  collocated 
DG  formulations  (where  the  interpolation  points  coincide  with  the  integration  points)  for  a 
specific  class  of  boundary  conditions  (see  [29]  and  [30]  for  the  solution  to  this  problem  for  the 
Navier-Stokes  equations). 

A  few  additional  comments  regarding  the  semi-implicit  discretization  are  in  order.  Since 
only  the  gravity  waves  (pressure  gradient)  are  discretized  implicitly,  with  the  Rossby  waves 
(advection  operator)  discretized  explicitly,  the  model  will  be  unconditionally  stable  with 
respect  to  the  gravity  waves  for  any  time-step  size  as  long  as  the  conditional  stability  with 
respect  to  the  Rossby  waves  is  satisfied  by  the  explicit  methods.  This  is  the  reason  why  we 
show  both  the  explicit  and  implicit  stability  regions  of  the  BDF  methods  in  Figs.  2a  and 
2b.  Of  course  one  could  also  choose  to  discretize  all  of  the  terms  implicitly  (including  the 
nonlinear  advection  operator)  -  this  is  known  as  the  fully-implicit  method.  The  reason  for 
choosing  the  semi-implicit  method  over  the  fully-implicit  method  is  that  in  doing  so  we  only 
need  to  contend  with  a  linear  matrix  problem;  for  the  fully-implicit  method,  we  would  have  to 
solve  a  nonlinear  matrix  problem,  which  while  not  impossible,  requires  many  more  iterations 
for  convergence  (outer  Newton  loops,  plus  inner  Krylov  loops,  [23]).  For  the  types  of  shallow 
water  problems  that  we  are  considering,  the  semi-implicit  method  should  be  more  efficient 
than  the  fully-implicit  method. 


5.  Numerical  Experiments 


For  the  numerical  experiments,  we  use  the  normalized  L 2  error  norm 


II  Ml’ 


ffi(0ex  act  0)  "dll 

inexact  dn 


computed  at  the  cubature  points  to  judge  the  accuracy  of  the  methods.  To  compute  the 
Courant  number  the  elements  are  decomposed  into  their  high-order  (HO)  grid  points  (which 
are  in  fact  the  Fekete  points)  and  these  grid  points  form  a  fine  grid  which  we  refer  to  as  the 
HO  cells.  The  velocities  and  grid  spacings  are  then  defined  at  the  centers  of  these  cells.  Using 
these  definitions  the  Courant  number  is  then  defined  as 

Courant  number  =  max  (  — — -  )  Ve  €  [1,  ...,Ne\  (30) 

V  As  /  HO 

where  C  =  U  +  \[§  is  the  characteristic  speed,  U  =  \Ju  ■  u  is  the  magnitude  of  the  velocity, 
and  As  =  y  Ax 2  +  Ay2  is  the  grid  spacing.  In  addition,  note  that  the  Courant  number  based 
on  the  advection  is  given  by  Eq.  (30)  with  C  =  U. 

We  use  the  symbol  nr  to  refer  to  the  refinement  level  of  the  grid.  This  variable  nr  represents 
the  number  of  quadrilateral  subdivisions  in  each  of  the  Cartesian  directions.  For  example, 
nr  =  1  corresponds  to  n'2r  quadrilaterals  and  2 n2  triangles;  the  factor  of  2  is  required  since 
each  quadrilateral  is  subdivided  into  2  triangles.  Examples  of  square  domains  with  nr  =  1, 
nr  =  2,  and  nr  =  4  are  shown  in  Fig.  3. 
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5.1.  Description  of  the  Test  Cases 

We  now  describe  the  test  cases  and  their  solutions.  It  should  be  noted  that  all  the  tests 
presented  below  require  no-flux  boundary  conditions  at  all  four  walls. 


5.1.1.  Linear  Standing  Wave  This  problem  involves  the  transient  solution  of  a  linear  inviscid 
standing  wave  without  rotation  which  sloshes  within  a  square  basin  of  unit  depth.  From  [22] 
we  take  the  analytic  solution  as 


h(x ,  t) 
u(x ,  t) 

v(x ,  t) 


cos  7 tx  cos  Try  cos  V^irt 
— —  sin  nx  cos  Try  sin  \Z2tt t 

V2 

— —  cos  ttx  sin  Try  sin  \z2nt 

V2 


with  (a;,  y)  £  [0,  l]2  and  t  £  [0, 2],  The  source  function,  S,  in  Eq.  (1)  is  zero  and  the  flux  tensor 
is  linearized. 


5.1.2.  Linear  Kelvin  Wave  This  problem  involves  the  transient  solution  of  the  linearized 
inviscid  equations  with  rotation.  From  [13]  we  use  the  analytic  solution 


h(x ,  t) 

=  1  +  exp 

u(x ,  t) 

=  exp 

v(x ,  t ) 

=  0 

2N  /  (a;  +  5-t)2 
exP - ^ - 


2N  /  (x  +  5-t)2 
exP - ^ - 


with  /o  =  0,  /3  =  1,  and  ( x,y )  £  [—10, 10]  x  [—5,5]  and  t  £  [0, 10]. 

5.1.3.  Nonlinear  Rossby  Soliton  Wave  This  problem  describes  an  equatorially  trapped 
Rossby  soliton  wave  [4].  The  soliton  wave  starts  off  in  the  center  of  the  domain.  It  then  moves 
westward  along  the  equator  without  changing  shape.  The  asymptotically  derived  analytic 
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solution  is  given  by 


h(x,y,t)  =  /i(0)+/i(1) 
u(x,y,t)  =  +u^ 

v(x,y,t)  =  v^+v^ 


where  the  superscripts  (0)  and  (1)  denote  the  zeroth  and  first  order  asymptotic  solutions  of 
the  shallow  water  equations,  respectively.  They  are  given  by 


h™  = 


u™  = 


.,(0)  = 


V 

dy 

di 


-9  +  6  y2 
4 

(2  y)e~^ 

3  +  6  y1 


and 


h =  c^rj^-  (— 5  +  2y2)  e  *  +  r?2$(1^(y) 
16 

=  c^r)-^-  (3  +  2y2)  e_J5"  -\-rf‘U^1\y) 
16 

u(1)  = 


where  y (£,f)  =  Asech2U^,  £  =  x  —  ct,  A  =  0.771  B2,  B  =  0.394,  and  c  =  c^0-*  +  c^1-*  where 
=  — y  and  c*'1^  =  —  0.395  B2.  The  variable  y  is  the  solution  to  the  equation 


dy 


any 


dy  d3y 


on 


Pn 


oa3 


=  0 


which  is  a  simplified  form  of  the  shallow  water  equations  upon  using  the  method  of  multiple 
scales  presented  in  [3].  Finally,  the  remaining  terms  are  given  by 


/  $(%) 
UW(y) 

V  Vd)(y) 


2  oo  /  (fin  \ 

=  e“T"  ^2  Hn(y) 

n— 0  \  Vn  ) 


where  Hn(y)  are  the  Hermite  polynomials  and  tpn,  un,  vn  are  the  Hermite  series  coefficients 
given  in  [4],  The  Coriolis  parameter  is  given  by  f(y)  =  y  where  (x,  y)  €  [—24,  +24]  x  [—8,  +8] 
t  €  [0, 10]  and  g  =  1. 

We  include  this  analytic  solution  for  completeness  but  one  cannot  use  this  test  for  validating 
the  exponential  convergence  of  the  methods  because  the  analytic  solution  is  only  a  first  order 
approximation.  However,  this  solution  can  be  used  to  check  the  phase  speed  of  the  soliton 
wave  simulated  by  the  numerical  model. 


5.1.4-  Linear  Stommel  Problem  The  linear  Stommel  problem  [34]  is  the  exact  steady-state 
solution  of  the  linearized  inviscid  equations  with  rotation,  wind  stress,  and  bottom  friction. 
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The  analytic  solution  of  this  problem  can  be  obtained  by  considering  the  linearized  momentum 
equation  as  follows 

V0  +  f(k  x  u)  +  7 u  =  r. 

If  we  define  the  Coriolis  term  as  f(y)  =  fo  +  (3  (y  —  and  the  wind  stress  as 


r  = 


i  +  0  j 


then  taking  the  curl  yields 


' )V2tp  +  (3 


dtp 

dx 


TqTT 

pHL 


where  we  have  written  the  velocity  field  in  terms  of  the  streamfunction  as  u  =  —  ^  and 

v  =  ■Ij.  Assuming  a  separation  of  variables  solution  of  the  type  tp(x,  y)  =  tp(x)  sin  yields 
the  following  second  order  ordinary  differential  equation  (ODE)  for  x 


TO 

pH' 


This  ODE  tells  us  that  we  need  to  seek  solutions  of  the  type  ip(x)  =  CeXx 
substituting  into  the  ODE,  gives  the  two  roots 


Cq  which,  after 


Al,2  — 


with  Co  =  -rjpff  ~  ■  Imposing  zero  streamfunction  boundary  conditions  at  the  domain  boundaries 
yields  the  final  solution  of  the  streamfunction  as 


where 


tP(x,  y)  =  (Co  +  CieXlX  +  C2ex*x)  sin  (^) 


Ci  =  C0 


1  -  eAai 


and 


C2  =  -C, 


1  —  e 


Ai  L 


°eA 2L  _  eA t_L  ’ 


This  then  yields  the  analytic  solution 


h(x) 

u(x) 

v(x) 


1 

9 


Co/3—  cos 

7 r 


(t)  +f(yW(x’y>> 


-pWcos(H) 

(C\\ieXlX  +  C2\2ex*x)  sin  . 


The  constants  required  to  completely  define  the  solution  are  /o  =  lx  10-4,  (3  =  lx  10— 11 , 
7  =  1  x  10-6,  g  =  10,  p  =  1000,  r  =  0.2,  H  =  1000,  and  L  =  lx  106.  The  models  are 
integrated  between  200  to  400  days  in  order  to  reach  steady-state.  We  regard  steady-state  as 
the  condition  when  the  error  norms  cease  to  decrease. 
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5.1.5.  Nonlinear  Riemann  Problem  We  used  the  nonlinear  Riemann  problem  previously  (see 
[20])  in  order  to  validate  the  spatial  operators  of  our  DG  model  and  its  slope  limiters.  We 
follow  the  outline  of  the  problem  presented  in  Toro  [37].  The  source  function  S  is  set  to  zero; 
this  leaves  a  balance  between  the  time  rate  of  change  of  the  conservation  variable  q  and  the 
the  divergence  of  the  flux  tensor.  Following  [37]  we  use 


h(x,y,  0) 


2.5  if  r<R 
0.5  if  r>  R 


with  u(x,y,  0)  =  0  V(x,  y)  €  [ — 20, 20]2  where  r  =  x 2  +  y2,  R  =  2.5,  and  t  €  [0,0.4],  The 
cylindrical  wave  is  positioned  initially  at  the  origin  and  moves  outward  for  increasing  time  t. 


5.2.  Comparison  of  the  Explicit  and  Semi- Implicit  Models 

In  the  sections  below,  we  compare  the  accuracy  and  efficiency  of  the  explicit  RK35  method  with 
the  semi-implicit  BDF  methods  of  order  I\  <  6.  For  all  simulations,  the  largest  Courant  number 
shown  for  RK35  represents  the  maximum  Courant  number  allowed  by  this  method.  The 
smallest  Courant  number  shown  for  the  semi-implicit  BDF  methods  represents  the  maximum 
Courant  number  allowed  by  the  explicit  BDF  methods;  the  only  exception  is  the  time-step 
convergence  study  that  we  now  describe. 


5.2.1.  Time-Step  Convergence  Study  The  first  study  we  conduct  is  the  rate  of  convergence  of 
explicit  and  implicit  time-integrators.  For  this  study  we  use  the  linear  standing  wave  problem 
with  nr  =  1  refinement  level,  corresponding  to  two  triangular  elements,  and  14th  order 
(. N  =  14)  polynomials.  For  this  resolution,  the  best  possible  normalized  hs  L2  error  norm 
that  can  be  achieved  by  the  model  is  1  x  10~9  which  we  obtained  experimentally  as  At  — >  0; 
this  we  consider  to  be  the  exact  numerical  solution. 

In  Fig.  4  we  show  the  results  for  the  RK35  and  BDF  methods.  The  maximum  time-step  used 
for  RK35  is  the  maximum  allowed  by  the  method.  For  the  BDF  methods  of  order  K  <  4,  the 
smallest  time-steps  corespond  to  the  explicit  methods  while  the  last  few  points  correspond  to 
the  implicit  methods.  For  the  BDF  methods  of  order  K  >  5,  all  the  simulations  are  obtained 
with  the  implicit  methods;  because  these  methods  are  high-order  in  time,  they  achieve  the  exact 
numerical  solution  for  time-steps  much  larger  than  those  allowed  by  the  explicit  method. 

We  define  the  time  rate  of  convergence  as 


rate  =  abs 


yl  log  [error Ati+1/errorAtJ  \ 

log  [Ati/Aii+i]  J 


where  At,  are  the  Nt  time-step  values.  Figure  4  shows  that  RK35  is  indeed  formally  third  order 
accurate,  and  that  the  BDF  methods  achieve  their  theoretical  values  of  order  K.  Furthermore, 
the  explicit  and  all  implicit  methods  yield  exact  mass  conservation  (to  within  machine  double 
precision)  not  just  for  this  test  case  but  for  all  cases  discussed  below.  Let  us  now  examine  the 
accuracy  and  efficiency  of  the  semi-implicit  BDF  time-integrators  for  various  test  cases  and 
compare  them  to  the  explicit  RK35  time-integrator. 


5.2.2.  Linear  Standing  Wave  In  Figs.  5a  and  5b  we  show  the  L 2  error  norms  and  the  wallclock 
time  as  functions  of  Courant  number  for  various  time-integrators  for  the  linear  standing  wave 
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Figure  4.  Linear  Standing  Wave.  The  normalized  hs  L2  error  as  a  function  of  time-step  for  various 
time-integrators.  All  runs  use  nT  =  1  and  N  =  14. 


a)  b) 

Figure  5.  Linear  Standing  Wave.  The  a)  normalized  hs  L2  error  and  b)  wallclock  time  as  functions  of 
Courant  number  for  various  time-integrators.  All  runs  use  nr  =  6  and  N  =  10. 


problem.  In  these  simulations  we  use  tenth  order  polynomials  (N=10)  with  72  triangular 
elements  (corresponding  to  nr  =  6).  Figure  5a  shows  that  the  RK35  explicit  method  is  more 
accurate  than  the  semi- implicit  BDF  methods  of  order  K  <  3.  For  I\  >  4  the  BDF  methods 
are  more  accurate.  Figure  5b  shows  that  all  of  the  BDF  methods  are  more  efficient  than  the 
explicit  RK35  method  for  the  same  Courant  number.  Furthermore,  the  BDF  methods  of  order 
K  <  6  allow  larger  Courant  numbers  than  the  explicit  RK35  method  and  the  lower  order  BDF 
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methods  are  more  efficient  than  the  high-order  BDF  methods. 

The  small  Courant  numbers  reported  for  BDF6  in  Fig.  5b  needs  to  be  explained.  To 
understand  these  results,  let  us  begin  by  discussing  the  stability  regions  of  the  implicit  BDF 
methods  in  Fig.  2b.  Along  the  imaginary  axis  all  the  BDF  methods  are  stable  for  large  z. 
Clearly,  for  some  range  of  z,  the  BDF  methods  of  orders  K  >  3  become  unstable.  This  is 
observed  in  Fig.  5b  for  Re(z )  =  0  and  \Im(z)\  <  5  for  BDF4,  for  example.  Once  this  first 
instability  region  is  reached,  we  stop  increasing  the  Courant  numbers  which  results  in  the 
small  Courant  numbers  reported  in  Figs.  5.  Note  that  we  do  this  for  all  of  the  simulations 
throughout  this  paper. 


5.2.3.  Linear  Kelvin  Wave  In  Figs.  6a  and  6b  we  show  the  normalized  hs  L2  error  norms  and 
the  wallclock  time  as  functions  of  Courant  number  for  various  time-integration  methods  for  the 
linear  Kelvin  wave  problem.  In  these  simulations  we,  once  again,  use  tenth  order  polynomials 
(N=10)  with  256  triangular  elements  (corresponding  to  =  16  and  n%  =  8).  Figure  6a  shows 
that  the  RK35  explicit  method  is  more  accurate  than  the  semi-implicit  BDF  methods  of  order 
K  <  3.  For  K  >  4  the  BDF  methods  are  more  accurate,  as  in  the  previous  case.  Figure  6b 
shows  that  all  of  the  semi-implicit  BDF  methods  are  more  efficient  than  the  explicit  RK35 
method  for  the  same  Courant  numbers.  Once  again,  the  BDF  methods  allow  larger  Courant 
numbers  with  BDF2  allowing  a  Courant  number  almost  one  order  of  magnitude  larger  than 
the  RK35  method. 


Courant  Number 


a) 


b) 


Figure  6.  Linear  Kelvin  Wave.  The  a)  normalized  hs  L2  error  and  b)  wallclock  time  as  functions  of 
Courant  number  for  various  time-integrators.  All  runs  use  n *  =  16,  ri%  =  8,  and  N  =  10. 


Because  the  BDF2  method  allows  such  large  Courant  numbers,  it  makes  it  very  difficult  to 
discern  the  performance  of  the  other  BDF  methods.  In  Figs.  7a  and  7b  we  show  the  results 
using  smaller  Courant  numbers.  In  these  figures  it  becomes  obvious  that  the  best  method  is 
the  BDF4  since  it  yields  the  best  error  norms  (together  with  BDF5  and  BDF6)  while  allowing 
for  larger  Courant  numbers  than  RK35,  BDF5,  and  BDF6  and  thereby  yielding  the  optimal 
result  when  considering  both  accuracy  and  efficiency. 
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a)  b) 

Figure  7.  Linear  Kelvin  Wave.  A  close-up  of  the  a)  normalized  hs  L2  error  and  b)  wallclock  time  as 
functions  of  Courant  number  for  various  time-integrators.  All  runs  use  n ^  =  16,  n %  =  8,  and  N  =  10. 


5.2-4-  Rossby  Soliton  Wave  Figure  8a  shows  the  wallclock  time  as  a  function  of  Courant 
number  for  various  time-integration  methods  for  the  nonlinear  Rossby  soliton  wave  problem. 
In  Fig.  8b  we  show  the  same  simulations  but  for  smaller  Courant  numbers.  Recall  that  for 
this  test  we  only  have  a  first-order  solution  which  is  not  sufficiently  accurate  for  performing  a 
convergence  study;  however,  we  can  use  it  to  determine  whether  the  solitons  are  moving  with 
the  proper  wave  speed.  We  consider  a  numerical  solution  to  be  accurate  if  it  agrees  exactly 
with  the  semi-analytic  solution  as  to  the  position  of  the  highest  peak  of  the  solitons. 

In  these  simulations  we  use  eighth  order  polynomials  (N=8)  with  384  triangular  elements 
(corresponding  to  n *  =  24  and  n-!  =8).  Figure  8  shows  that  the  semi- implicit  BDF  methods 
K  <  6  are  more  efficient  than  the  explicit  RK35  method  for  the  same  Courant  numbers. 
Additionally,  the  semi-implicit  BDF  methods  K  <  3  admit  larger  Courant  numbers  than  the 
explicit  RK35  method.  Since  this  case  is  nonlinear,  both  the  explicit  and  implicit  BDF  methods 
are  used  in  tandem  to  solve  the  problem.  Therefore  in  this  test  case,  the  stability  regions  of 
both  the  implicit  and  explicit  BDF  methods  are  relevant.  Looking  at  the  stability  region  of  the 
explicit  BDF  methods  given  in  Fig.  2a  we  note  that  the  BDF5  and  BDF6  have  particularly 
small  stability  regions  that  result  in  the  small  Courant  numbers  reported  in  Fig.  8  for  these 
methods.  The  BDF  methods  of  order  K  <  4  have  larger  stability  regions  in  both  the  explicit 
and  implicit  forms  and  is  the  reason  why  these  methods  perform  more  efficiently  than  the 
high-order  BDF  (K  >  5)  and  RK35  methods. 

5.2.5.  Linear  Stommel  Problem  Figure  9a  shows  the  wallclock  time  as  a  function  of  Courant 
number  for  various  time-integration  methods  for  the  linear  Stommel  problem;  in  Fig.  9b  we 
show  a  close-up  of  the  same  simulations  for  smaller  Courant  numbers.  For  this  case  we  have 
a  steady-state  analytic  solution  and  so  the  accuracy  of  the  time-integrator  only  plays  a  small 
role.  The  accuracy  of  the  model  is  completely  dependent  on  the  polynomial  order  of  the  DG 
method;  the  only  role  that  the  time-integrator  has  is  to  maintain  the  stability  of  the  solution 
while  doing  so  as  efRcienctly  as  possible. 

In  these  simulations  we  use  eighth  order  polynomials  (N=8)  with  32  triangular  elements 
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a)  b) 

Figure  8.  Rossby  Soliton  Wave.  The  wallclock  time  as  a  function  of  Courant  number  for  various  time- 
integrators.  Figure  a)  shows  the  results  for  large  Courant  numbers  while  b)  shows  them  for  small 
Courant  numbers.  All  runs  use  n *  =  24,  n%  =  8,  and  N  =  8. 


(corresponding  to  nr  =  4).  Figure  9a  shows  that  the  BDF  methods  K  <  2  allow  larger 
Courant  numbers  than  the  other  methods.  More  importantly,  Fig.  8b  shows  that  all  of  the 
BDF  methods  are  more  efficient  than  the  explicit  RK35  method  for  the  same  Courant  numbers. 
In  addition,  for  this  test  case,  the  implicit  BDF  methods  admit  as  large  a  Courant  number  as 
the  explicit  RK35  method.  The  fact  that  the  implicit  BDF  methods  are  more  efficient  than 
the  explicit  RK35  even  for  the  same  Courant  number  is  impressive  especially  since  the  implicit 
BDF  methods  require  much  more  machinery  to  solve  the  problem.  Recall  that  implicit /semi- 
implicit  methods  require  the  use  of  iterative  solvers  (in  this  case  GMRES)  and  preconditioners 
(in  this  case  Jacobi  preconditioning)  in  order  to  solve  the  resulting  linear  matrix  problem.  Even 
with  all  of  this  machinery,  the  implicit  methods  are  more  efficient  than  an  explicit  RK  method 
-  this  does  not  seem  possible  at  first  glance.  The  reason  for  these  surprising  results  is  simple: 
for  Courant  numbers  less  than  1  the  implicit  BDF  methods  require  fewer  than  5  GMRES 
iterations  to  converge;  recall  that  the  RK35  method  requires  5  stages.  Thus  at  this  range  of 
Courant  numbers  the  implicit  BDF  methods  are  more  efficient  than  RK35  with  respect  to 
operation  count  which  translates  to  smaller  wallclock  times.  For  the  larger  Courant  number 
values,  the  number  of  iterations  are  greater  than  5  but  the  larger  time-steps  (hence  fewer  time- 
integration  loops)  compensate  for  the  extra  costs  incurred  with  respect  to  operation  count. 

One  final  comment  is  in  order.  Since  this  test  case  is  linear  and  the  solution  is  steady-state, 
we  could  have  used  infinitely  large  Courant  numbers  for  the  implicit  BDF  methods  K  <  2. 
We  have  only  chosen  to  report  the  maximum  Courant  numbers  that  maintained  stability  for 
the  nonlinear  Stommel  problem.  Since  we  do  not  have  an  analytic  solution  to  the  nonlinear 
Stommel  problem,  we  then  use  the  linear  problem  to  ensure  that  we  are  achieving  L2  errors 
of  1  x  10~6  which  is  the  exact  numerical  solution  for  eighth  order  polynomials  with  nr  =  4 
(see  [20]).  This  means  that  the  results  reported  in  Fig.  9  are  also  representative  of  the  types  of 
efficiency  gains  offered  by  the  semi-implicit  BDF  methods  for  nonlinear  problems.  However,  it 
should  be  noted  that  the  reason  why  such  large  Courant  numbers  can  be  used  in  this  test  case 
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a)  b) 

Figure  9.  Linear  Stommel  Problem.  The  wallclock  time  as  a  function  of  Courant  number  for  various 
time-integrators.  Figure  a)  shows  the  results  for  large  Courant  numbers  while  b)  shows  them  for  small 
Courant  numbers.  All  runs  use  nr  =  4  and  N  =  8. 


has  to  do  with  the  disparity  between  the  speed  of  the  gravity  waves  (the  height  of  the  fluid) 
and  the  Rossby  waves.  For  the  Stommel  problem,  the  gravity  waves  are  much  faster  than  the 
Rossby  waves  and  this  will  be  the  case  for  all  deep  ocean  flows.  Let  us  now  discuss  a  type  of 
flow  problem  for  which  the  semi-implicit  method  is  not  well  suited. 

5.2.6.  Nonlinear  Riemann  Problem  In  Figs.  10a  and  10b  we  show  the  wallclock  time  as  a 
function  of  Courant  number  for  various  time-integrators  for  the  nonlinear  Riemann  problem. 
Note  that  we  only  report  the  explicit  RK35  and  explicit  BDF  methods.  The  results  for  this  test 
show  that  the  explicit  BDF  methods  of  order  K  <  4  compete  with  RK35  in  terms  of  efficiency. 
The  stability  regions  of  the  BDF  methods  of  order  I\  >  5  are  too  small  and,  while  faster 
than  RK35  for  a  given  Courant  number,  cannot  compete  with  the  maximum  Courant  number 
admitted  by  RK35.  Let  us  now  discuss  why  we  do  not  show  results  for  the  semi- implicit  BDF 
methods. 

We  cannot  use  the  semi-implicit  BDF  methods  for  this  case  because  the  Rossby  waves  are 
faster  than  the  gravity  waves.  This  means  that  the  linearization  used  to  construct  the  equations 
in  (6)  is  no  longer  valid.  The  linearization  used  in  the  current  semi-implicit  formulation  assumes 
that  (f>B  is  much  greater  than  <ps  which  is  not  true  for  the  Riemann  problem  (as  is  evident  by 
the  initial  conditions  where  (j>s  =  0.5  and  max{(j>s)  =  2  meters2 /second).  We  show  the  result 
of  the  Riemann  problem  only  to  point  out  the  limitation  of  our  current  approach.  Let  us  now 
discuss  some  possible  solutions  to  this  dilemma. 

Defining  the  Froude  number  as  the  ratio  of  propagation  speeds  of  Rossby  (call  them  R)  and 
gravity  (call  them  G )  waves,  then  if  the  flow  is  subcritical  (i.e. ,  G  >  R)  then  the  fix  to  the 
problem  is  relatively  simple.  Instead  of  linearizing  about  a  constant  state,  say  4>b  ,  we  linearize 
instead  about  the  known  state  at  the  current  time-step  (i.e.,  +  <^g,  where  n  denotes  the 

current  time-level).  This  presents  very  little  changes  to  the  current  semi- implicit  approach. 
On  the  other  hand,  if  the  flow  is  supercritical  (i.e.,  R  >  G)  then  nothing  in  the  semi-implicit 
machinery  can  improve  the  efficiency  since  the  terms  responsible  for  the  fastest  waves  in  the 
system  are  discretized  explicitly  in  time. 
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The  simplest  solution,  given  the  methodology  described  in  this  paper,  is  to  switch  from  the 
semi-implicit  to  the  explicit  methods  which  is  achieved  by  setting  the  parameter  Ssi  =  0  in 
the  code  -  this,  of  course,  has  to  be  done  with  the  additional  constraint  that  the  time-step  be 
changed  in  order  to  satisfy  the  explicit  stability  region  of  the  BDF  methods. 


a)  b) 

Figure  10.  Riemann  Problem.  The  wallclock  time  as  a  function  of  Courant  number  for  various  time- 
integrators.  Figure  a)  shows  the  results  for  large  Courant  numbers  while  b)  shows  them  for  small 
Courant  numbers.  All  runs  use  nr  =  100  and  N  =  1. 


Another  approach  is  to  discretize  the  equations  fully  implicitly  in  time  which  then  requires 
the  solution  of  a  nonlinear  matrix  problem.  We  prefer  the  first  approach  for  its  simplicity 
and  to  this  end  we  are  developing  tools  to  automate  the  selection  of  the  time-step  as  well 
as  the  value  of  the  switch  Ssi ■  Our  analysis  has  shown  that  not  all  K  order  methods  are 
created  equally.  For  example,  taking  all  the  results  collectively  shows  that  the  RK35  method 
behaves  most  like  the  BDF4  method  and  so  the  optimal  combination  would  be  to  use  the  semi- 
implicit  BDF4  as  long  as  the  linearization  is  valid  and  then  switching  to  the  explicit  RK35 
when  the  linearization  breaks  down  or  supercritical  flow  is  encountered.  The  value  of  such  a 
hybrid  solution  strategy  can  be  appreciated  by  considering  the  semi-implicit  time-integration 
of  a  tsunami  wave  beginning  in  the  middle  of  the  deep  ocean.  As  the  wave  approaches  the 
coastline,  the  semi-implicit  linearization  breaks  down  and  the  flow  becomes  supercritical  which 
then  requires  the  code  to  switch  to  explicit  mode.  We  hope  to  report  the  results  of  such 
simulations  in  the  near  future. 


6.  Conclusions 

We  present  a  high-order  family  of  semi-implicit  time-integration  methods  based  on  backward 
difference  formulas  (BDF)  for  the  triangular  discontinuous  Galerkin  method  as  applied  to  the 
oceanic  shallow  water  equations;  We  use  a  high-order  discontinous  Galerkin  method  defined 
on  unstructured  triangular  elements  which  is  especially  useful  when  attempting  to  resolve  the 
complex  geometry  resulting  from  the  representation  of  coastlines  in  coastal  ocean  models.  In 
this  work,  we  have  extended  the  explicit  in  time  high-order  DG  method  that  was  shown  to 
be  exponentially  convergent  (for  smooth  problems)  to  semi-implicit  in  time.  The  semi-implicit 
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BDF  time-integrators  of  order  K  >  4  are  shown  to  yield  better  accuracy  than  the  third  order 
explicit  Runge-Kutta  method.  Furthermore,  the  BDF  methods  of  order  K  <  4  require  far 
less  wallclock  time  to  deliver  these  solutions.  We  show  that  the  semi-implicit  BDF  methods, 
even  without  any  optimization  and  without  the  use  of  sophisticated  preconditioners,  yields 
better  efficiency  than  the  explicit  RK  method.  We  expect  that  on  a  parallel  computer,  with 
the  aid  of  preconditioning  and  reduction  of  the  implicit  problem  to  a  Helmholtz  problem,  the 
speed-up  of  the  semi-implicit  method  compared  to  the  fastest  explicit  methods  will  be  further 
increased.  In  future  work  we  plan  on  adding  lateral  diffusion,  variable  bathymetry,  and  wetting 
and  drying  algorithms  to  the  model  in  order  to  perform  tsunami,  storm  surge,  and  inundation 
simulations. 
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